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Abstract 

At asymptotically high energies, thermalization in heavy ion collisions can be described via 
weak-coupling QCD. We present a complete treatment of how thermalization proceeds, at the 
parametric weak-coupling level. We show that plasma instabilities dominate the dynamics, from 
immediately after the collision until well after the plasma becomes nearly in equilibrium. Initially 
they drive the system close to isotropy, but Bjorken expansion and increasing diluteness makes the 

-12 , 

system again become more anisotropic. At time r ~ ce~^~Qj the dynamics become dominated by 
a nearly-thermal bath; and at time t ~ a~Q~ l the bath comes to dominate the energy density, 
completing thermalization. After this time there is a nearly isotropic and thermal Quark-Gluon 
Plasma. 



I. INTRODUCTION AND SUMMARY 



Relativistic heavy ion experiments have given us new insights into the dynamics of hot 
QCD [l], |2|]. One of the most striking results is the success of a hydrodynamic description 
[l, 0]. This seems to imply that the hot QCD matter generated in a heavy ion collision 
thermalizes quite quickly and then has a period as a nearly thermal medium well described 
by hydrodynamics (which is predicated on the assumption of local near-equilibrium). Theo- 
retically, we have no idea why this should be true. In fact, we do not even understand early 
thermalization in the theoretically clean case where we consider the limit of large nuclei at 
arbitrarily high energy per nucleon. In this limit, the very early dynamics are believed to 
be well understood, in the framework of the "colored glass condensate" 0, |6[ and its post- 
collision debris, the "glasma" 0. According to this description, the initial conditions should 
be approximately boost-invariant, intense classical colored fields with a single energy scale 
Q s , the "saturation scale". Boost invariance means that the natural coordinates are the 
transverse directions x, y, rapidity rj, and proper time since the collision, r. We will rename 
t to t and generally work in terms of z = rjr, the physical distance along the beam direction 
in T), t coordinates, which expands (red-shifts) linearly with time (Bjorken expansion). The 
scale Q s increases with increasing nuclear size and with increasing energy per nucleon 0, 
so in principle this scale becomes large in the limit of interest (though whether this is true 
in practice at achievable energies is another matter). Provided that Q s 3> Aqcd, physics at 
the scale Q s should be characterized by weakly-coupled QCD, with the coupling taking a 
characteristic value a s (Q s ) <C 1 (henceforth we write a s (Q s ) as a). Weak coupling should 
make the theoretical problem easier to address. Nevertheless, the subsequent evolution of 
these fields is not well understood. 

The most comprehensive attempt to date to describe thermalization at central rapidities 
from these initial conditions was by Baier, Mueller, Schiff and Son (BMSS) They argued 
that the medium made up of p ~ Q s excitations becomes anisotropic and dilute; but LPM 
suppressed radiation of lower energy daughters generates a bath of p <C Q s excitations 
which come to dominate the dynamics and lead to the breakup of the p ~ Q s "hard" 
"parent" excitations into a thermal bath. However, as pointed out by Arnold, Lenaghan 
and Moore their discussion ignored the physics of plasma instabilities [HI, T2], which 
should dominate at least the early stages of the dynamics. The role of plasma instabilities in 
these early stages was considered analytically by Bodeker [l3| and by a mix of analytic and 
numerical considerations by Arnold and Moore • Classical field simulations with plasma 
instabilities were considered by Romatschke and Venugopalan IE)], and several groups made 
related studies [l6|. However, none of these studies followed the dynamics past the early 
stages or explained how final thermalization occurs. 

In this paper we will take advantage of our recent study of the dynamics of weakly 
coupled gauge theories out of equilibrium (Ref. (T3], henceforth "KM1") to give a purely 
analytical and parametric description of how "the glasma" should thermalize when Q s is 
such a large scale that « <C 1. The initial state of the glasma should be well described by 
intense classical fields with much longer coherence in the r/ direction than in the transverse xy 
plane [3] • That is, the initial conditions are very anisotropic, with a longer coherence length 
(smaller mean wave number) in the 77 direction. But we find (as sugg ested by Ref. 0) that 
plasma instabilities rapidly reorganize this into a configuration which is only anisotropic at 
the 0(1) level. 

From this point we find two attractor solutions. In one solution, the plasma becomes 
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highly anisotropic. Bjorken expansion increases the anisotropy but does not lower the wave- 
number of typical excitations. This leads to a dilute plasma of wave-number p ~ Q s excita- 
tions. This solution is similar to the one found by BMSS 0, except that plasma instabilities 
dominate the dynamics (as proposed in Ref. 10] and discussed in Refs. [H, 14]). In this 



solution the momentum space anisotropy grows with time as £s. In the other attractor solu- 
tion, the plasma becomes nearly isotropic at a time t longer than Q~ x by at most logarithms 
of a. Plasma instabilities resist Bjorken expansion to keep the plasma nearly isotropic. The 
typical excitation energy falls with time, and the typical occupancy falls but not as fast as in 
the first solution. This solution is somewhat like the one recently suggested by Blaizot, Gelis, 



Liao, McLerran, and Venugopalan |18j, except that we find that plasma instabilities play a 
key role in the dynamics and that no Bose-Einstein condensate forms. Fully nonperturbative 
and nonequilibrium early-time dynamics determine which attractor the system finds. But 
we argue that the path from the initial conditions to the more anisotropic attractor is much 
more likely. 

Assuming the first (strongly anisotropic) attractor is correct, we find that a nearly- 
thermal bath of low-momentum (p ~ T <C Q s ) excitations forms and comes to dominate 
the dynamics at time t ~ a~^~Q~ 1 . By t ~ a^Q~ x the soft bath causes the "hard" 

p ~ Q s excitations to disintegrate and the system becomes nearly isotropic and thermal, 

—11 — —i 

with a temperature Tfi na i ~ a^rQst~ ■ But while nearly thermal, the system is initially 
non-Newtonian; anomalous viscosity [l9[ (due to plasma instabilities) controls the level of 
isotropy until t ~ a~Q~ x , after which the fluid becomes Newtonian. Before t ~ a~Q~ x 
the plasma is far from equilibrium; we parametrically estimate the (quite complicated) phase 
space particle number density at all momenta p < Q s and all times from Q~ l to final ther- 
malization. The results are presented in Figure |5] and Table [III 

Our results indicate a faster breakup of the hard p ~ Q s excitations than the old estimate 
t ~ a^Qj 1 from BMSS |||, but still somewhat slower than the absolute lower bound of 



a a Q s 1 derived by Arnold and Lenaghan [20|. But it is consistent with the estimate for 
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full thermalization times presented in KM1, who claim t cq ~ a~QJT~; inserting Tg na i 



returns t ~ a 2 Q~ s 



II. INITIAL CONDITIONS AND ATTRACTORS 

As we briefly mentioned in the introduction, an extensive literature on the small- a; behav- 
ior of large nuclei (see Ref. [1]) and small-a; evolution in heavy ion collisions (!, 0] indicates 
that, at time t ~ Q s , the debris of a heavy ion collision is dominated by classical gluon 
fields with coherence length I ~ 1/Q S in the xy plane, much longer coherence in z = rjt, and 
energy density e(t ~ Qs 1 ) ~ • At times t > Q~ l these fields have lost phase coherence 

and the evolution should be described with the tools we introduced in KM1, modified to 
include the Bjorken expansion of the geometry. 



A. Combining KM1 with Bjorken expansion 

We begin by introducing notation, including the notation used in KM1. In KM1 we 
considered the evolution of a system with a single dominant characteristic wave number 
Q in a time-independent geometry. The system is characterized by a scale Q (which need 
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Variable 


Describes 


Detailed Definition 


a 


time 


/ „ — a/O— 1 


b 


Typical momentum 


Q ~ a b Q s , where f(p > Q) « /(p ~ Q) 


c 


Mean occupancy 


/(p = Qx) ~ a~ c 


d 


Anisotropy 


d < : f(Qx) - f(Qz) ~ a~ d /(Q) 

d>0: f(p z ,p±) <C /(0,p ± ) i f > a d |P±l 



TABLE I: Brief summary of our log notation for the most important descriptors of the system 



not be Q s ] we will write Q = a b Q s ), a typical occupancy f(p ~ Q) ~ ct~ c , and a degree 
of anisotropy characterized by d. For d > most excitations have p z < a d p and we write 
5 = a d . The occupancy f(p) ~ a _c applies for p z < Sp, outside this range the occupancy is 
small, so the particle number and energy densities are 

for d > 0, n ~ c/" c Q 3 and e ~ a d " c g 4 ~ a d - c+4b Q A s . (2.1) 

For d < the particle number is nearly direction independent, but the direction dependence 
varies with relative amplitude e = a~ d - that is, (f(Qx) — f(Qz))/f(Q) ~ e. Therefore, for 
d < 0, the particle number and energy densities are 

for d < 0, n ~ a~ c Q 3 and e ~ a" c g 4 ~ a~ c+ib Qt . (2.2) 

We will also introduce one more logarithmic variable; we will write the time as t = a~ a Q~ 1 
(or a = \n(tQ s ) / hi(l / a)) . The initial conditions are fixed at t ~ which is a = 0. 

Parametrically later times are a > 0. Our notations are summarized in Table [H 

According to KM1, the occupancy-anisotropy or c, d plane is divided into four regions, 
shown in Figure [TJ At large anisotropy the dynamics are dominated by magnetic or "Weibel" 
~2l| plasma instabilities. That means that certain long-wavelength magnetic fields grow to 



large amplitudes; there are associated electric fields but they are much weaker. Magnetic 
fields change an excitation's direction without changing the energy; so the main physics in 
this region involves the randomization of excitations' directions. Induced (bremsstrahlung) 
emission or merging of excitations also occurs but is subdominant. For less isotropic but 
over-occupied systems (1 > c > and d < (c — l)/3), ordinary scattering dominates the 
dynamics. Elastic scattering drives the system towards isotropy, and inelastic scattering 
drives the mean occupancy c towards zero. The two types of scattering are of comparable 
efficiency (in stark contrast to scalar field theory). For underoccupied systems below some 
level of anisotropy, the main physics is that a bath of smaller-momentum excitations forms 
and eventually causes the dominant excitations to shower and join this bath. The last region, 
of extremely high occupancies, exhibits Nielsen-Olesen instabilities 22[. This region has no 



relevance for us. The arrows in Figure [TJ show qualitatively how these dynamics cause c, d 
to evolve with time. 

All this is for a homogeneous and non-expanding system. But the debris of a heavy ion 
collision expand. In our choice of coordinates the system is nearly statistically homogeneous, 
but the z direction grows linearly with time. This Bjorken expansion red-shifts away the 
z component of momentum, increasing the anisotropy of the system, and raising d. And 
for d < it lowers the typical excitation's energy, increasing b. In the scattering, plasma 
unstable, and Nielsen-Olesen regions the evolution of b, c, d will be determined by an interplay 
of the equilibration processes and Bjorken expansion. 
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FIG. 1: Simplified version of Figure 1 from KM1 [l7j |. The vertical axis represents how anisotropic 
the dominant excitations are; the horizontal axis shows how large the dominant excitation occu- 
pancies are. The labels indicate what the dominant physics is in each region; the arrows suggest 
how those dynamics "move" the system in this two-dimensional space. There are no arrows in 
the "soft particle bath" region because here the typical excitations break up and are replaced by 
something very different. 



B. Initial conditions and attractors 



The initial conditions we described above have d > and c = 1 + d at time t ~ Qj 1 - This 
is the point marked "Initial Condition" in Figure [2j The figure also shows two attractor 
solutions which we find, which we now explain in more detail. By an attractor we mean 
a path on the c, d plane, such that a system relatively close to the path will fall onto and 
then follow the path as a result of the combination of dynamics and Bjorken expansion. 
There are two attractors because there is a fundamental difference between the dynamics of 
a system with large anisotropy versus a system with small anisotropy. This difference turns 
out to allow solutions of either type. 

Consider first a system which is highly anisotropic, so the typical excitation has p z <C 
Then Bjorken expansion does not change the typical excitation's energy. And angle 
randomization - such as that induced by plasma instabilities - increases the region of phase 
space which is occupied, and therefore brings down the typical occupancy (reduces c). But 
consider instead a system which is only weakly anisotropic. Typical excitations have \p z \ ~ 
\p\, and therefore Bjorken expansion reduces typical excitation energy. And since angle 
randomization is mixing modes of nearly the same occupancy, it has little impact on the 
typical occupancy. These fundamental differences mean that we have to treat the two cases, 
d > and d < 0, separately. And we will find an attractor in each. Specifically, we will find 
for d > 0, 

a T cl 

d——, c=l , b = 0, (attractor 1) (2.3) 
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FIG. 2: Occupancy-anisotropy plane, the initial conditions from heavy ion collisions, and the two 
attractor solutions we find. 



while for d < there is another attractive solution, 

56a , 31a 



d 



1 - 



(attractor 2) 



135 ' 135 ' 135 

These attractors, and the initial conditions, are displayed in Figure [2j 



(2.4) 



C. First (anisotropic) attractor 

First consider the case d > 0, where most excitations have p z <C \p\. This region lies 
entirely in the plasma instabilities dominated region in Figure [1] provided c > —1/3 (as it 
initially is). To describe evolution in this region it is convenient to think of the descriptors of 
a system, 6, c, and d, as functions of time or of a, and to ask about their a derivatives; d = 
dc(a)/da and similarly for b', d'. Each quantity evolves due to Bjorken expansion and due to 
the effects of dynamics, principally plasma instabilities; we will write d = c'| B j orken + c'| dyn 
to distinguish the part of the evolution from each effect. 

The effect of Bjorken expansion is a red-shifting of the p z component of momentum, 



^jj^- = — ^ (Bjorken expansion) . 



(2.5) 



This does not change c; occupied states retain their occupancy, they just change the value of 
momentum. It also does not change 6; since p z <C \p\, the typical excitation's energy is not 
changed. Thus, c'| B j orkcn = = &'| B j rken- Bjorken expansion does lead to an increase in 
d. Recall that 5 = (p z /p) RMS = oi d . The angle 8 evolves according to 

do" do" dp z 5 dp z 5 
dt dp z dt p z dt t ' 

which using Eq. ( 12. 5 j) and the definitions of a, d becomes 

d(d) d(d)dt d(ln(5)/ln(a)) dt 



(2.6) 



d! 



tdS 



Bjorken 



da 



dt da 



dt 



d(ln(tg s )/ln(a- 1 )) 5dt 



+1. (2.7) 



6 



Now we need to add the effect of plasma instabilities. As mentioned, the dominant thing 
the instabilities do is to deflect excitations, widening their angular range. This changes c 
and d, but not the typical energy b; b'\ dyn = 0. Since &'| B j orken = as well, we find b' = 
and therefore b = 0; the typical excitation energy remains Q = Q s . Dynamics also cannot 
change the total energy density e. According to Eq. (12. ip . e depends on the combination 
d — c + 4b. Since b' = 0, we learn that a"| dyn — c 'ld yn = 0- Combining with Eq. ( 12. 7p . we 
learn that d' — d — 1. Since initially c — d = 1, we have 

b — 0, c - d = 1 - a (d> case) . (2.8) 

Therefore the state of the system at a given time t (or a given value of a) is determined by 
d(a); c is fixed to be c = 1 + d — a. 

To finish determining the dynamics we need to find out how fast plasma instabilities cause 
the level of anisotropy to change. According to KM1, there is a time scale tbroadcn which 
determines how long it takes the plasma instabilities to significantly increase the range of 
angles of p ~ Q excitations (to reduce d). The anisotropy therefore evolves at 

— — => d'l = = -a- a Q:H- 1 j . (2.9) 

A+ + 'dyn ± ^ s broaden \ / 

Clf- ''broaden ^broaden 



According to KM1 jlTj] , the time scale ^broaden is 



^broaden ~ { 3c+£,dS ^ , , ' (2.10) 

a 2 Q 1 d < ^ . 
Combining with Eq. ( 12. 8 P and Eq. ( 12. 9p . we find 

At a time a during the evolution, d might take its value on, above, or below the attractor. 
Since we have two constraints, Eq. ( I2.8p . on b, c, d, the possible deviation from the attractor 
is completely fixed by the value d(a). If d(a) > a/8 then d'\ dyn <ti —1 and d falls very 
rapidly, on a time scale shorter than the age of the system. If d < a/8 then d'\ dyn ~ and 
dl = 1. In this case d grows linearly with a, rapidly exceeding a/8. Therefore the solution 
d = a/8 is an attractor. The rest of Eq. (I2.3P follows from Eq. (12. 8p . According to Figure HI 
this solution makes sense until it reaches the point d = 3/10, c = —11/10, which occurs at 
a = 12/5 (t ~ a ~Q^ x )- After this time a soft bath forms and takes over the dynamics. 
We will return to these "late-time" dynamics in the next section. 



D. Second attractor 

Now consider the case d < 0, that is, systems which are nearly isotropic. Now Eq. ( 12. 5ft 
is still correct, but since (p 2 z ) ~ p 2 /3, expansion now changes the typical particle energy, 

^lBjor k en = l/3. (2.12) 
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Correspondingly, the energy density falls with time as e ~ t 4 / 3 (t 1 due to dilution and 
t -1 / 3 due to the work done by longitudinal pressure). Then Eq. ( 12. 2 p becomes (compare 
with Eq. (EE}) 

-c' + 46' = 4/3 or 4b = c-l + Aa/3. (2.13) 

We can use this expression to eliminate b in favor of c, d in what follows. 

Next we find the effect of expansion on d. A reduction of p z by p z — > p z (l — e) would 
induce an anisotropy of e. Since d = ln(e)/ ln(l/a), d is very sensitive to expansion: 

^Bjorken ~ ^ = > 1 . (2.14) 

In any self-consistent solution, d'\ dyn will have to be compensatingly large, so d! = (d'\ B ^ orkcn + 
d'\ dyn ) £ I- 

We can immediately rule out the existence of an attractor solution in the scattering region, 
d < (—1 + c)/3. In this region, number-changing processes are as efficient as direction- 
changing elastic processes. 1 Therefore c'| dyn ~ ^1d y n- 2 So in the scattering region, either 
d'\ A ~ 1, in which case d rises very fast due to d'\ B ^ olken ^> 1, or d'\ dyn ^> 1, in which 
case c'| d ^> 1 and the system evolves very fast to lower occupancy. Since lower occupancy 
features slower dynamics, both c'| dyn and d'\ dyn would then collapse, again allowing d to rise 
into the plasma-unstable region. 

So consider the region with important plasma instabilities, > d > (— l + c)/3. According 
to KM1, there are two relevant time scales here: tbmaden, the time scale for plasma instabil- 
ities to cause large-angle direction change; and i mer ge, the time scale for plasma- instability 
induced number- changing processes. (Think of these as hard synchrotron radiation, in- 
duced as the excitations bend in the magnetic fields; but for over-occupied systems the 
synchrotron absorption process is more important than emission.) The time scales are (also 
using Eq. (|2TT3|) ) 

C — 1 3d-3+3c _-, 3rf-3+3c-26+2q 18d+15(e-l) + 8a 

> d > : tbroadcn ~«2Q~tt 2 t ~ a 12 t, 

3d-7+7e i 3d-7+7e-4b+4a 9d+ 18(e- l) + 8a 

Emerge ~a 4 Q ~« * t~OJ 12 t. (2.15) 

Merging (two particles joining into one with the sum of the initial particles' energy) raises 
the typical particle energy and lowers the occupancy. So 

0>d> C —±: c >\ dyn ~--^. (2.16) 

<J " merge 



1 Since there is some confusion in the literature on this point 18(, we will summarize why this is true. The 
rate for an excitation to undergo a generic-angle 2 <H- 2 process is a 2_2c Q _1 (a 2 for the two vertices, 
oT 2c for two Bose factors). Small-angle scattering is enhanced by a factor Q 2 /m 2 ~ a~ 1+c because of 
the Coulomb divergence; but a small angle scattering does little to drive the system towards equilibrium. 
Adding a hard radiation or merging onto a small-angle scattering process reduces the rate by a factor of 
a 1 ~ c (a for the vertex, a~ c for a Bose stimulation factor for the new external state; phase space factors 
are canceled by collinear factors as usual in ISR or FSR), yielding a rate of a 2_2c Q _1 . If f(p <C Q) oc p" 1 
then the rate in the infrared is actually larger by a factor of Q/p, which makes the development of a 



negative chemical potential impossible (contrary to the assumptions of la]). This point is discussed at 



2 



some length in Ref. [23[, section HIE. 

The only exception is if c = 0, so number changing processes are in equilibrium. That is, there can be 
nearly-thermal, scattering dominated systems, just not ones with paramctrically large occupancy. 



s 




FIG. 3: Evolution near the second attractor. At a given time a, one point on the attractor line 
(circled in blue) is preferred. Above the solid blue line, the dynamics rapidly lowers the anisotropy; 
above the line, Bjorken expansion rapidly raises it (as suggested by green arrows). On the blue 
line, number-changing processes move the system towards the blue circled point. 



Broadening leads towards isotropy, lowering d. But so does merging; as discussed in 
KM1, when merging is induced by plasma instabilities, the particles near the xy plane 
merge more often than those with large p z , reducing the number of in-plane excitations. 
This last effect is independent of the anisotropy e, while angle randomization only affects 



the O(e) angle- anisotropic part of the distribution. Therefore [17 



<*v=-(V^+*vM • (2 - i7) 

\ ''broaden ''merge / 

The tbroadcn term dominates if > d > (— 1 + c)/7; for d < ( — 1 + c)/7 the t mer ge term 
dominates. 

Now we show that Eq. (12. 4p constitutes an attractor. While in the d > case we had 
two constraints on b, c, d, now we have only one constraint, Eq. (12.131) . So if we are off the 
attractor at time a, we are at some point in the c, d plane. Figure [3] shows how c, d might 
deviate from the attractor solution, Eq. (12. 4p . at some given time a. The blue line is the 
line where -d'\ dyn = d'| Bjorken : using Eq. (12141) . Eq. fl^TTjl and Eq. fl2~T5l) . 

f d — — ^ d > 

- d'\ - d'\ if I 2 15' 7 ' (2 181 

d Idyn " d iBjorken \ = 2(1 - C ) - f d < . 1 j 

Above this line, |d'| d | > d'\ B ^ olken . Since both are large, d! is large and negative here. 
Below the line, |ef| d | < d'\^ or]ken and d! is large and positive. In either case the dynamics 
drive d to the line on a time scale short compared to the system's age. 

Next we must check the evolution of c on the blue line in Figure [3j to see whether it 
is driven towards the attractor. To the left on the line, where d > <i a ttract and c < extract, 
Eq. (I2.15p . Eq. (I2.16p . and Eq. (I2.18P show that d ~ 0. Since the attractor moves towards 
lower c with time (Eq. (I2.4p ). points here move towards being on the attractor. To the right 
of the line, where c > extract, Eq. (I2.16P and Eq. ( I2.17P show that c'| dyn ~ Q^'ldyn ~ — 1- 
In fact it is more negative than -1, since every merging process changes typical occupancy, 
but only the anisotropic part of merging processes change anisotropy. Therefore c would fall 
quickly when it is greater than the attractor value. This verifies that the solution, Eq. ((27 
constitutes an attractor. 
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FIG. 4: Zoom-in on the region of the c, d plane near c = 1, d = 0, showing the early-time evolution 
from the initial conditions. It takes some time for the instabilities to develop, during which Bjorken 
expansion increases anisotropy. When plasma instabilities finally act, they bring the system to a 
region where the two attractors are well separated, leading to attractor 1. 



E. Why we think the first attractor is realized 

The initial conditions are c = 1 + d and d > 0, corresponding to an energy density 
s ~ a~ 1 Q^. Using Eq. ( j2.1ip . as soon as a > this initial condition should fall rapidly to 
approximately the point c = 1, d = 0. That is the point where the two attractors meet. So 
which attractor will the system actually follow? 

Here we argue that, if we try to include corrections suppressed by a logarithm of a, it 
appears that the first (anisotropic, d > 0) attractor is preferred. Consider first the initial 
conditions. The characteristic wave- number (or inverse coherence length) is Q s , the energy 
density is e ~ a~ 1 Ql/t. Under these conditions plasma instabilities grow at a characteristic 



rate called m, the screening scale [12|; for these conditions m 2 ~ ae/Q 2 s ~ Q s /t. The 
unstable modes have p z 3> p±, a region which does not receive large occupancy in the initial 
conditions. Instead, these plasma-unstable modes presumably start with typical occupancy 



/ 



2 ■ 



the amount provided by vacuum fluctuations. The plasma-unstable fields have little 



impact on the dynamics until they have grown to occupancy ~ a . This requires a time 
of order t ~ Qj 1 ln 2 (l/a), or a ~ 21n(ln(l/a))/ln(l/o;). Therefore, there is a short delay 
before the plasma instabilities start to do their job. During this time, Bjorken expansion 
raises d, moving it away from the d = c — 1 line. So when the system starts to become more 
isotropic, it does so with a mean occupancy below c = 1, arriving away from the c = 1, d — 
point. We illustrate this idea in Figure HJ By the time d approaches 0, the two attractor 
lines are separated, and we naturally land from "above" (the anisotropic side) on attractor 
1. 



3 Some readers might worry that it is cheating to seed a classical field with initial conditions set by the 
size of quantum fluctuations. For a field which experiences an exponential instability, it is not cheating, 



a point discussed at some length in the inflationary preheating literature, see for instance Ref. 24 
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Of course, it might somehow occur that the amplitude of plasma instabilities overshoots 
the level where nonlinear interactions normally limit instability growth. That would allow 
the unstable fields to grow larger, and angle change to be faster, than the estimates given 
above (which are based on steady-state plasma instability evolution }17|). For instance, it is 
known that, if the system is very aniso trop ic and if instabilities start with extremely small 



initial amplitudes, this can occur |25l . |26[. But we claim that even such an "overshoot" 
cannot lead to near-isotropy. 

First, note that the plasma- unstable fields have much larger magnetic than electric field 
amplitude. Therefore, they deflect hard excitations more efficiently and they absorb their 
energy. The time scale for large-angle deflection is shorter than the time scale for large energy 
absorption; so the plasma-unstable fields never dominate the energy density, which remains 
in the hard excitations. So the only way to get isotropy is to get large-angle deflection. 

We showed above that the scales m and Q s become separated; m <^ Q s (at least by 
log(l/ct)). The unstable fields have wave numbers k ~ m. For a field of wave number 
k ~ m <C Q s to bend excitations of momentum Q s by and order-1 angle in a time scale 
1/m, the force would have to be dp/dt ~ mQ s , requiring a magnetic field of strength 
B ~ mQs/ g. But in nonabelian field theory, whenever a magnetic field of wave number k 



exceeds the field strength B > k 2 /g, it becomes Nielsen-Olesen unstable 22|, leading to the 



very rapid (t < 1/k = 1/m) collapse of any such magnetic field. Therefore the magnetic 
fields associated with plasma instabilities can never grow larger than B ~ m 2 /g - too small 
to deflect hard p ~ Q s excitations by large angles. 

In the remainder of the paper we will concentrate on filling in all details about the 
development of the system, assuming it follows the first attractor. 



III. DETAILED EVOLUTION ALONG THE FIRST ATTRACTOR 

In the previous Section, we described how the occupancies c and the anisotropy parameter 
d of the primary particles with p ~ Q s evolve with time along the attractor 1, until at the 
time a = 12/5, when the momentum diffusion experienced by the primary particles is taken 
over by the soft sector. In this Section we describe the details how the soft sector evolves 
during this time, and we describe in detail what happens after the soft sector takes over 
the dynamics. In particular, we find the occupancies and angular ranges of all excitations 
throughout the evolution of the system. The results of this Section are collected in Figure 
Sand Table M 

We divide the evolution in three main stages. In the early stage the primary particles have 
occupancies larger than 1, and all processes involving these modes are stimulated. During 
this stage, the system grows steadily more anisotropic and radiated daughters become less 
and less important. Incidentally, this means that classical Yang-Mills theory under Bjorken 
expansion never isotropizes and the interactions become weaker with time. 

At the time a = 8/7, the attractor 1 crosses the c-axis in Figure [21 and the occupancies 
of the hard excitations become small. Therefore, the Bose stimulation of processes involv- 
ing hard excitations is lost. This is the middle stage. The loss of stimulation makes the 
interaction rates decrease more slowly, and the soft sector will start to affect the dynamics. 
The soft sector is also anisotropic and therefore has its own unstable modes. During the 
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FIG. 5: Log-log plot of the time-momentum plane, divided into regions with distinct behaviors for 
the angular range of excitations about the xy plane, 9, and the occupancy f(p). The horizontal 
axis is a defined as t = a~ a Qj 1 or a = — ln(iQ s )/ln(a). The vertical axis is the log of momentum. 
The equations describing lines are written in shorthand, so for instance (8 — 3a)/16 means p = 
q sQ! (8-3o)/16 The i etters A-T refer to Table [Til 

second stage, the momentum diffusion 4 q = dp\/dt arising from these new plasma instabil- 
ities increases relative to the "primary" instabilities of the hard particles. The new plasma 
instabilities come to dominate q at large angles at the time a = 96/55. But q at narrow 
angles, experienced by the primary excitations, remains dominated by the "primary" insta- 
bilities. Most of the particle number in the soft sector is in a nearly isotropic bath, whose 
distribution forms an f{p) (xl/p "tail" below a cutoff scale p max - Initially f(p ma , x ) S> 1, but 
at the time a = 56/25, f(p mayL ) = 1 and the soft bath becomes nearly thermal. 

Plasma instabilities from the nearly-thermal bath grow to dominate at all angles at time 
a = 12/5, which is where attractor 1 in Figure |5] intersects with the boundary between 
the regions labeled "Plasma Instabilities" and "Soft Particle Bath Forms". Thereafter, the 
broadening of the hard primary particles will be controlled by the soft sector and its plasma 
instabilities. The instabilities cause hard particles to undergo successive splitting processes, 
depositing their energy into the thermal bath, increasing its temperature. At time a = 5/2, 
the primary p ~ Q s particles have had time to break up completely, leaving only a nearly 
isotropic thermal bath. We now discuss these stages in more detail. 



4 Here and in the following, "transverse" or p± refers to the directions transverse to the initial momentum 
of the particle, not the directions transverse to the beam axis. 
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TABLE II: Angular range of excitations and typical occupancies for the regions labeled in Figure 

El 

A. Early stage: a < 8/7 

1. Angular distribution, the scale k\ so 

We will now find the angular ranges and occupancies of the soft sector during the evo- 
lution. Plasma instabilities give rise to angle dependent transverse momentum diffusion 
characterized by q of size 17 



q 



5 2 m 3 ~ a°* Ql 



Pz ^ $P± , P > 5 2 m ~ oaQ s , 



5-i#-i m 3 ^ B~ x a—Ql , p z /p ± = 9, p> m/ {89) 
where m is the thermal mass 



m 2 ~ 



d 3 p 



« / —f(p) ~ c^-^Q ~ a a Q 
p 



m ~ a a/2 Q s . 



The angular range which excitations of momentum p occupy will be given by 
9V ~ q(9)t ~ onrB- 1 ar a Q 2 s 6 ~ a*(Q a /pfi . 

We define the scale where 6* ~ 1 as fci SO : 

1 ~ at(g s //c iso )l =4> fc iso ~ a^Q s . 



(3.1) 



(3.2) 



(3.3) 



(3.4) 



Below this scale particles are nearly isotropic. 



2. Particle production 



The instabilities not only change the angular range the soft particles populate; they 
also affect their occupancies through instability induced splitting processes. Provided p > 
5~ 2 m ~ a^Qs, particle emission is LPM suppressed. The formation time for an emission is 



ti^(p,9)^q(e)p-^a^Qlp-\8l9) 



(3.5) 
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but for our case, in a time tf orm the radiated daughter only accumulates a transverse mo- 
mentum of V qt{ orm < Sp, so 9 ~ 5. Ignoring stimulation, the rate each hard particle splits 
off daughters of momentum p is dT/d\n(p) = r sp i it ~ G^fem- Summing over the number 
density hard particles nhard and including the final-state hard particle stimulation factor 
[1+/(Q S )] ~ oT c ^ the number of daughters radiated per logarithmic momentum range in 
time t ~ a~ a <57 1 is 

ridau g htcr(p, a) rsj n haxd a~ c at^J: ~ a d - c Q 3 s a 1 ~ c a^'a~ a {Q s /p)^ ~ a~ 1+ ^Q]p~^ . (3.6) 
3. Redshifting and old particles 

3 

As a consequence of expansion, the rate of particle production falls quickly (oc a^ a ), and 
therefore the distribution of soft particles at a given time a need not to be dominated by 
the particles emitted during that time scale. The daughters radiated at an earlier time a' 
are still around, though they have reduced in number and energy due to Bjorken expansion. 

For p > fc iso , the longitudinal momentum component is parametrically small compared 
to the total momentum, and the effect of the expansion is simply to reduce the number of 
particles by a factor a a ~ a . So the number density of particles, emitted at time a' , which are 
present at time a > a' is 

^redshifttjo > &iso, a; a') ~ q"""' n daughtcr (p, a'). (3.7) 

For p < k lso , the expansion also decreases the z-component of momentum. The plasma 
instabilities keep the soft distribution below /c; so isotropic so that the net effect of the ex- 
pansion and angle randomization is to redshift the isotropic spectrum below A; iso to smaller 
momenta, much like in the case of attractor 2; 

n T edsUft(p < kso, a; a') ~ « a_a ' ^daughter («"^ (a " a 'V, o!) . (3.8) 

Both above and below /c iso , the particle number density at momentum scale p is dominated 
by particles which where emitted at the scale k lso . For p > k lso , the dominant contribution 
to n(p) arises from time a' when fci SO was p, namely a a ~ {p/Qs)~ , so that the particle 
number density above A; iso is 

n(p) ~ a °- a V 1+ § a 'Qfp-5 rsj a- 1+a (p/Q s )fQ 3 sJ (3.9) 
and the corresponding occupancies in the occupied angle region are 

f(p)~^~*- 1+la (Qs/p) l «- (3-10) 

This is Region A in Figure EJ 

For p < fciso, the particles created at the scale fci SO at time a' have just redshifted to scale 
p at time a provided that 

a^- a ')fc iso (a')~p, (3.H) 
so that the particle number density at time a is dominated by particles created at the time 

a a ' ~ a f a (Q s /p)f. (3.12) 



14 



The particle number density and occupancies then read 



n(p) ~ a^'a-^'QUa-^-^p)-^ a- 1+ ^ a (Q s /p)^Ql (3.13) 



n(p) 



27, 



f(p) ~ ~ a- 1+ T3 a (g s /p)i3. (3.14) 



This is Region B in Figure El 



4- Saturation, the scale p max 

As seen in KM1, distributions that have infrared tails steeper than that of a thermal 
distribution rearrange themselves very effectively to form a thermal-like tail so that 

f(p) oc l/p, for p < p max . (3.15) 

In particular, if the occupancies of daughters f(p) exceed 

f{p) ~ (3.16) 

the rate for the daughter particles to rejoin the hard particles becomes larger than the 
corresponding emission rate, and the soft particle number density in fact decreases due to 
inelastic scatterings. This is the case in the Region B below p ~ a^Q s . The total rate for 
soft particles rejoining to hard particles is (see Eq. (2.54) of KM1) 

r rej oin(p) ~ 4t / d " k T sMp)lf(k)f(p) ~ nY T d Tsplit 2 {p) ~ a^Q s {Q s /p)i , 

n{p) J k ~Q s k Q s p 2 

(3.17) 

so that the modes at scales below p max , with r rc j in(Pmax)^ ^ 1 or 

Pmax ~ (3.18) 

have had time to change their occupancies by an order 1 amount and adjust to the form of 
Eq. (I3.16p . Below p max the occupancies then read 

f(p)~a- 1+7 s%Q s /p). (3.19) 

This is the region C of Figure [51 



B. Middle stage: 8/7 < a < 12/5 
1. a > 8/7 New particles, the scale k Q \d 

For a > 8/7, the Bose stimulation for the particle production is lost. As a result, the 
particle production rate decreases more slowly and particles produced at later times may 
start to dominate the particle number density. The number of soft particles created at the 
time scale a is now (remove ct~ c from Eq. (13. 6p ) 

^dau g htcr(p, a) ~ n Ylwd atJ^ m t ~ afQ 3 s (Q s /p)^ . (3.20) 
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For k > k- lso this becomes dominant over the redshifted particles for k < /c Q id given by 
(compare Eq. ([3T20]) with Eq. <K9\f ) 

Km ~ a*-*Q s . (3.21) 

For k < k lso the new particles become more numerous for k > k D \d with (Eq. (13 .201) and 
Eq. flSH) 

Kid ~ a~^ +7 ^Q s . (3.22) 
These are the regions D and E in Figure [51 respectively. Region B ends at a = |f. 



2. a> 16/11: Soft joining, the scale p max 

We defined the scale p max as the one below which f(p< p max ) oc 1/p. At early times we 
saw p max is set by the rate of emission/absorption off hard p ~ Q s excitations. At later times, 
merging processes between p ~ p max excitations can become faster and come to control this 
scale. The rate for soft joining (in Region E, just above p max ) is 

r mergc (p) ~ ai^fo 6 ~ l)[l+/(p)] ~ « 1+ ^(Q S /#Q, , (3.23) 

and r merge (p max )t > 1, using p max from Eq. (13.181) . after time a = 16/11. After this, the 
evolution of p max is controlled by soft merging; 

r mcrgc (p max )t ~ 1 => p max ~ a3+Mg s . (3.24) 

Modes with p < p max effectively thermalize with each other. The energy density is domi- 
nated by p ~ p max modes, so the occupancy is determined by the energy density of particles 
residing at the scale p max . Hence, for p < p max , 

Sip) « 8 128 Q s /P - 3.25 

p 

This is Region G in Figure [5j 



3. a > 96/55: iVew; plasma instabilities, the scale k new .i ns t 

The distribution of daughter particles is anisotropic; therefore the daughter particles also 
give rise to plasma instabilities. As discussed in KM1, the strongest instabilities are driven 
by daughters at the scale fc iso , where the angular distribution is order-1 anisotropic. The 
resulting plasma instabilities give rise to momentum diffusion 

<?new-inst ~ ^{kso) ~ Qft/t (k iso ) kf so , (3.26) 

which is comparable in all angular directions 9 (see Eq. (4.21-4.22) of KM1). By comparing 
Eq. ( 13. II) to Eq. ( 13.261) . using Eq. ( 13 .41) and Eq. (13.201) . we find that g ne w-inst = <?inst at large 
angles 9 ~ 1 at a = 96/55. After this, the evolution of the scale k iso is driven by the 
instabilities caused by the scale h lso itself. The scale ki SO can be found by self-consistently 
solving 

k lo ~ <?ncw-inst£ and g ncw _ ins t ~ f^(k iso )kf so with /(fc is0 ) ~ a^Qi^J (3.27) 
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resulting in 

24-a 24+33a , -42+23a 

kso ~ a 68 <2s , <?ncw-inst ~ « 34 , /(«iso) ~ « 34 ■ (3-28) 

When g ncw -inst dominates, particles at scales above /c; so will have an angular range 9 ~ 
a~^~ (Qs/p) ■ The range of angles where g new -mst > <7inst is # > ex we . Using Eq. (13.31) . 
we see that excitations experience g nC w-inst provided p < A; ncw _ inst ~ a ^~^Q S . This is the 
region H in Figure 

The new q changes the rate of the soft merging and therefore the functional form of p max 
as well as f(p < p max )- Repeating the analysis of Subsection II II B 21 with q ~ g nC w-inst gives 

184+15a 605q-920 

Pmax ~ a 5« Q s an d f(p) ~ a ioss (Q s /p). (3.29) 
This is the region I in Figure [5j 



4- a > 56/25: Thermalization of the soft sector and the scale k sp m 

At the time a = 56/25, f(p max ) drops to ~ 1 and the cascade of the infrared modes 
by joining processes, which started at the time scale a = 16/11, is complete. The soft 
sector now forms a nearly-thermal bath with a temperature T ~ £> m ax- Subsequently, the 
dominant inelastic process changes from joining to splitting, and those particles with p > T 
that have had time to undergo a democratic splitting (i.e., emit a daughter with a momentum 
comparable to emitter's momentum) in less than the system age cascade down in energy by 
multiple splittings, and lose their energy to the thermal bath. The modes which have had 
time to undergo a democratic splitting are those for which p < k sp \[ t , with 

. . I r> o 92 — 35a . . 

r sp lit(fcsplit)t ~ O^form^ ~ 1 =^ ksplit ~ Ot ^ncw-inst^ ~ Ot 34 Q g . (3.30) 

The temperature of the bath is set by the highest momentum modes falling onto it, so that 
the temperature becomes 

T ~ e\ ~ (/(A; split )A; s 4 plit )3 ~ oT^Q s . (3.31) 

This explains Region J in Figure |5j 

Most of the excitations below ki SO but above T have had time to split and release their 
energy to the thermal bath. However, a small portion of the particles, the ones created 
in the last time scale ^residence ~ r split (p) ~ (p/& S plit) ^> are still around and have not 
yet undergone the cascade to lower scales. The flux of energy from the scale fc sp iit moving 
through scales to the scale T is 

f t ~ e(k split )/t ~ f(k spht )kt p Jt, (3.32) 

so that the energy density carried by particles created in the last t res idence time scale at the 
momentum scale p is 

de 

e {p) ~ Residence "77, (3.33) 

at 

so that f(p) and n(p) for /c iso > p > T are 

f(p) ~ e(p)/p' ~ f(k spm )(k sp]it /p)l ~ a*{Q s /p)l, n(p) ~ a"(Q s /p)*Q* (3.34) 
which is a continuation of Region E. 
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5. a > 160/69: Instabilities from the nearly thermal bath 

At time a = 160/69 the scale k sp nt crosses the scale k [so . This reduces the number of 
excitations at the scale fci so , which ceases to dominate as a source of plasma instabilities. 
However, the thermal bath is incompletely isotropic, so it also gives rise to plasma insta- 
bilities, with associated momentum broadening parameter gVinst- Beyond a = 160/69, gYinst 
dominates for 9 ~ 1. 

Two mechanisms keep the soft nearly-thermal bath incompletely isotropic. First, Bjorken 
expansion continually "stretches" the thermal bath. Secondly, because g ne w-inst is order-1 
anisotropic, the radiated daughters which arrive at the scale T do so with an 0(1) small-# 
bias. Both mechanisms return a level of anisotropy for the bath of 

T 2 

Cbath ~ 7 , (3-35) 

gTinst£ 

which is the ratio of a typical p\ in the bath to the change in p\ during the system's age. 
Using Eq. (3.16) of KM1, we find 

<?Ti„st ~ 4ath^(^) 3 ~ ef ath a§T 3 • (3.36) 
Substituting gVinst into Eq. f)3.35p . we find 

3 2 3 12 3 

ebath ~ «"s(tT)"5 , g Tinst ~ asT^fs . (3.37) 

The temperature is still determined by the k sp nt scale, T ~ (/(^ S piit)) 1 ^ 4 ^spiit, and the splitting 
scale is given by fc sp ii t ~ a 2 grinst£ 2 - Solving self-consistently for e, gnnst, k Bp nt, and T gives 

13__3o_ A 12 i 15q a 26_41a -44+23a 

T~«2s ii2Q S) qrinst^a 7 2S Q S , fepiit ~ a 7 28 Qs , fbath ~ a ™ . 

(3.38) 

This explains Region K in Figure |5j 

The new dominant q changes the rate of angular broadening and hence also the scale 
^new-inst- The momentum range where g ne w-inst dominates, and the angular range of excita- 
tions, are 

36-15a 48-13a 

grinst > ginst if P < ^new-inst ~ a ™ Q s where 9 ~ a se Q s /p. (3.39) 

This explains Regions L and M, above and below k^, respectively. 

Below /c S piit but above T the cascading particles start their fall in a small angle range 
at the scale fc sp ut and as they cascade down in energy, their angular range broadens. The 
cascading particles stay at the momentum scale p for the time ^residence ~ (p/^spiit) 1 ^ 2 ^ and 
during this time their angular range grows to 

p 2 2 ~ gTinst^residence 9 ~ tt^^ (Qs/p) * (3.40) 

with 9 ~ 1 for p < k- lso ~ a 8 s4 J Q s . The estimate for the particle number density from 
Eq. (I3.34p still holds so that the occupation number below fc split is f(p) ~ (Qs/p)^ ■ 

This is Region N in Figure |5j 5 



5 A quick calculation of q due to instabilities associated with the scale fci so returns <7 n ew-inst ~ a 58 Qs ~ 
gainst- All scales from fci so down to T actually have comparable roles in setting q; but this does not change 
parametric estimates. 
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C. Late stages: 12/5 < a < 5/2 



At a = 12/5 we find gi nst ~ (ZTinst- Beyond this scale, the new plasma instabilities from 
the incompletely isotropic thermal bath dominate momentum broadening at all scales, i.e., 
this is the first moment when the soft sector starts to influence the dynamics of the hard 
particles. This is the point d = 3/10, c = —11/10 in Figured where attractor 1 enters the 
region marked "Soft Particle Bath Forms." As a result, the anisotropy of the hard modes 
departs from d = a/8 to become more isotropic, and the rate of radiation of daughters 
increases, since it is now controlled by the larger grinst- 

The splitting scale is now given by fc sp iit ~ « 2 ^rinst^ 2 - Since there is only a single q, k sp u t 
is the scale where there is of order 1 daughter per hard particle. So 

T 4 ~ rihard/csput ~ a 1_a gTinst<5 s (3.41) 
while Eq. ( 13 . 3 7|) still holds. Solving solving self-consistently, we find 

grinst ~ a 3 Q 3 s , T ~ a 1_ tQ s , fc split ~ a 5 ~ 2a Q s , e bath ~ . (3.42) 



This is Region O in Figure [5j 

The new q broadens the angular distribution more effectively, and for p > k sp \[ t we 
find 6(p) ~ a^Qs/p. In particular, the angular range of the primary hard excitations is 
5 ~ a~z~ , which grows with time for the first time since a ~ 0. However, the hard excitations 
remain anisotropic to the bitter end. 

Repeating the calculation of the behavior of particles within the cascade, we get the 
results for regions P, Q, R, and S in Figure [5j 

The thermal bath continues to draw energy from the hard excitations until k sp n t reaches 
Q s and the hard excitations, which dominate the system's energy, themselves split and join 
the thermal bath. This occurs at a = 5/2. At this stage the temperature is T ~ a^Q s , and 
the residual anisotropy of the nearly-thermal bath is e ~ a*. This is still larger than the 
value e ~ as where a thermal bath becomes dominated by scattering rather than plasma 
instabilities. 

A shortcut calculation which leads to the correct final equilibration time is the following. 
The energy density scales with time as e ~ a a ~ 1 Q^, so the would-be equilibration temper- 
ature is T(a) ~ a^Qg. According to KM1 17j], the equilibration time for an anisotropic 
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system with T < Q s is 

t eq ~ a^rQlT^r . (3.43) 

(This result was derived by asking how long it takes a thermal bath at temperature T to 
cause hard splitting in excitations of momentum Q s , given that the thermal bath will be 
somewhat anisotropic and the dynamics are dominated by plasma instabilities; see KM1 
[13].) Equating t eq with the age of the system t ~ a~ a Q~ x and inserting our expression for 
T, we find a = 5/2 or t eq ~ a^Qj 1 . 



D. Aftermath: a > 5/2 

After a = 5/2 the p ~ Q s initial excitations are gone, leaving behind a nearly thermal 
bath. The system is again characterized by a single scale, the temperature T, and hence 
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the system can again be described using Figure [2j At the time a = 5/2, d — — 1/4 and 
c = 0; subsequent evolution makes the distribution more isotropic, and the solution will 
move down along the <i-axis. 

The system is now nearly isotropic, and therefore the energy density begins to scale with 
time as e oc t _4//3 rather than as t~ l . Using the value of the energy density at a = 5/2, we 
find 

T~af-Mg s . (3.44) 

Initially q remains dominated by plasma instabilities from the residual anisotropy of the ther- 
mal bath giving rise to "anomalous viscosity" . We still find Eq. ( I3.37P to hold; substituting 
the known T value, 

7a _1_ o 4a 5_ / \ 

q ~ a 5 2 Q S , e ~ a 15 12 . (3.45) 

This drops below the value e ~ a 1 / 3 , entering the regime where ordinary scattering dom- 
inates, at a = 45/16. Beyond this time scale, the system behaves like an ordinary ther- 
mal bath expanding with a viscosity given by the normal (non-anomalous) viscous law; 
q ~ a 2 T 3 ~ a a+ iQ z s and the residual anisotropy is e ~ a^~^. 

In comparison, attractor 2 reaches the rf-axis at the time a = 135/56, so that the sys- 
tem becomes nearly thermal slightly earlier than in the case of the first attractor. As a 
consequence of being nearly isotropic from the beginning, the temperature is now only 

T at t 2 ~ «MQ s (3.46) 

resulting in 

7a o 4a 1_ . . 

9att2 ~ « 5 Qs > 6 att2 ~ « 15 2 • (3.47) 

In the case of attractor 2, scattering comes to dominate at the time a = 25/8, after which 

<?att2 ~ a a+ 'iQl and e att2 ~a^". 

In both cases, the behavior after scattering comes to dominate continues, in principle, 
indefinitely. At exponentially late times, a ~ 3a _1 //?o ln(a _1 ), the temperature scale reaches 
Aq CD and hadronization occurs. 

IV. DISCUSSION 

We have presented a complete parametric description of thermalization after a heavy 
ion collision in the theoretically clean, in-principle limit of arbitrarily large and high-energy 
ions, such that the scale Q s (which dominates particle production from the original collision) 
is taken so large that a s (Q s ) 1. Our most significant finding is that thermalization 
does eventually occur. Plasma instabilities dominate the dynamics at all times until well 
after the system becomes nearly isotropic and thermal. Initially, they cause the highly 
anisotropic starting configuration to become order-1 isotropic. But Bjorken expansion, and 
falling typical occupancies, allow the system to again become anisotropic - though the level 
of anisotropy only grows as t». 

The process of thermalization is dominated by 

1. the generation of plasma instabilities due to the anisotropy of the typical p ~ Q s 
excitations; 

2. the radiation of daughter excitations, which eventually form a nearly-thermal bath; 
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3. the plasma instabilities generated by the residual anisotropy of these daughter ex- 
citations, which are responsible for finally breaking up the p ~ Q s excitations and 
thermalizing the system at a time scale t ~ a~^Q~ x . 



After near-thermalization the system remains dominated by plasma instabilities of a nearly- 
thermal bath (anomalous viscosity) until the time scale t ~ a~^Q~ l , after which elastic 
scattering becomes the most important physical process and the plasma is well described by 
standard kinetic theory 

Since there is an extended period of expansion in which the equation of state is e oc 
t~ rather than £ oc t"J as for a thermal system, there is significant entropy generation. 
Most of the entropy of the final system is generated around t ~ a~iQ~ l , when the initial 
hard excitations fragment and join the thermal bath. This entropy generation could have 
phenomenological consequences in connecting the Colored Glass initial conditions to final 



state multiplicities, see [27 



It is not clear to us how high Q s would actually need to be for the parametric estimates of 
this paper really to make sense. However we feel that this work does resolve the long-standing 
puzzle of what happens in a heavy ion collision in what should be the most theoretically 
clean limit - that of extremely large saturation scale. 

It would be nice to replace this parametric estimate with a numerical estimate, by com- 
puting the "constant" in tthcrm = The good news is that this only requires one 
to study the final stage of thermalization, 12/5 < a < 5/2 (since all other stages take para- 
metrically less time). To study the relevant physics in this stage, we would need to determine 
the relationship between anisotropy e and momentum diffusion q for a nearly thermal bath, 
and we would need to compute the exact rate of fragmentation for a hard excitation experi- 
encing a given value of q. The former could be done by studying the plasma instabilities of 
a weakly anisotropic system. This should be doable using existing methods which combine 
classical field dynamics with anisotropic hard loops. The latter is also straightforward. 
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